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Abstract 

The motion of three-dimensional (3D) solitary waves and vortices in nonlinear crystal-like struc- 
tures, such as photonic materials, is studied. It is demonstrated that collective excitations in these 
systems can be tailored to move in particular directions of the 3D system. The effect of modula- 
tion instability is studied showing that in some cases it can be delayed by using a lensing factor. 
Analytical results supported by numerical simulations are presented. 

PACS numbers: 05.45.Yv,42.65.Tg,03.75.Lm 
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I. INTRODUCTION 



Recently, great progress has been made in the experimental fabrication of three-dimensional 
(3D) crystal-like optical nano structures, such as photonic crystals 1] and photonic meta- 
materials j^]. Light in 3D crystal-like optical nano structures can interact with the regular 
pattern of the structure setting up resonances. These resonances can cause beams and 
pulses to be deflected in unconventional directions and even to slow down the speed of light. 
Engineered optical metamaterials with unique electromagnetic properties, have become in 
recent years a hot research topic because of their interesting physics and exciting potential 



applications as, e.g cloaking |3j, |4(, or negative refractive index |5[], among others. 

The theoretical analysis of these systems have been mostly performed with the help 
of numerical simulations, where the 3D spatial distribution of the effective electromagnetic 
properties of the material medium (such as the permittivity and permeability) are tailored to 
have specific electromagnetic properties in the continuum limit 3j. Other discrete effects of 
the system on the propagation of light waves are usually neglected. Similar approximations 
have been also adopted for studying cloaking of matter waves 0]. So far, most of the studies 
have been done for linear planewaves, so nonlinear effects have been neglected. There is, 
however, strong evidence fro,n .ow-dnnens.ona. systen. flflflQ that d.screte effects in 
combination with nonlinearities may play an important role in describing realistic 3D crystal- 
like structures. So, nonlinear excitations such as moving solitary waves (in short solitons) 
and vortices can be expected. 

With respect to the theoretical models, it has been shown that electromagnetic waves 
interacting with metamaterials in the tight binding limit can be described as a 3D lattice of 
microresonators modeled by the 3D discrete nonlinear Schrodinger equation (3D-DNLSE) 
j^]. In the case of 3D photonic crystals, so far, no discrete model has been proposed. 
However, it is well known that the dynamics of light beams in photonic waveguides can 
be described by the 2D-DNLSE jsl, 3|, where the discreteness is transversal to the light 
propagation. Moreover, it has been also shown that trapped light waves travelling along 
chains of optical resonators can be described by the ID DNLSE [y, |7|, where the discreteness 
is longitudinal to the light propagation. These two low-dimensional models strongly suggest 
that 3D photonic crystals with high index of refraction may be effectively modeled 

in the tight binding limit by a 3D lattice of optical resonators governed by the 3D-DNLSE 
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FIG. 1. (Color online) Superposition of snapshots of p m ,n,p isosurfaces at different time values of 
a soliton-soliton-soliton collision. The isosurfaces are defined as Kp max (t), where the maximum 
Pmax{t) = max mi „ jP Pm,n,p{t), and k = 0.5. The arrows (in yellow) are meant to guide the eye and 
show the path and direction of motion of the solitons. At t = [labels A, B, and C (red color)] 
the initial conditions follow from Eq. (JSJ), where (A) k = 0.95§{0, 0, -1}, r m>Tl)P = {32,32,52}, 
n = 1.05, (B) k = 0.95f {-1, -1,0}, r mi „ iP = {52, 52, 0}, n Q = 1.1, and (C) k = 0.95f {1, 1, 1}, 



L m,n,p 



{12, 12, 12}, f^o = 1.5. At t = 10 [label D (green color)] the collision can be observed. 



And at t = 24 [labels A', B', and C (blue color)] the solitons after the collision can be observed. 



Other parameters are U = —1, J x 
so only sech-type solutions in Eq. 



J y = Jz = 1, lx = ly = Iz = 1/3. Here, (3 X = (3 y = Z = +1, 
\ are considered. 
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121. 



Notice that the 3D-DNLSE is an ubiquitous dynamical-lattice model which may emerge 
from a variety of other important problems and has its direct physical realization in Bose- 
Einstein condensates (BECs) trapped in strong optical lattices jlOl. I 111. Il2|. 

The aim of the present work is to study the dynamics of moving solitons and vortices in 
nonlinear 3D-crystal-like structures described by the 3D-DNLSE. 
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FIG. 2. (Color online) a: Evolution of the maximum p max {t) = max mjn ^ Pmnp(f) f° r t ne collision 
in Fig. [D (Black solid line). For comparison, p ma x(t) of the individual solitons (cases where the 
solitons do not collide) plotted in Fig. [1] is presented (red, blue, and green dashed lines; lines lie 
very near each other), b: some solutions of the Eq. (|20p for the cases (3 = 1 [u = 0, (black solid 
line), v = 1 (red dashed line), v = 2 (blue dotted line and green dot-dash line)], and (3 = —1 [v = 1 
(orange double-dot-dash line)]. 

II. THEORY 



The general form of the 3D-DNLSE with coupling constants J x , J y and J z is 

m,n,p—l 4~ 1pm,n,p-{-l) 
U Pm,n,p 'lprn,n,p 0, (1) 

where 

Pm,n,p = IV'wijTijpl (2) 

can be interpreted as a intensity (e.g. in crystals built of microresonators) or as a probability- 
density (in BEC arrays). In Eq. ([1]) the nonlinear coefficient U is a real constant and t is 
the time coordinate. 
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In order to proceed we shall consider an expansion into harmonics of a travelling wave 
ansatz for an envelope complex function, i.e 

CO 

1pm,n,p = ^ X^( S AW?) eX P( i6 W)> ( 3 ) 
/x,i/,£=l 

where 

= r - v^t (4) 

and 

= k /*.",f ' r ~ ( 5 ) 

Here, r = {m,n,p} is the position vector, v M)l/ ^ = {i^, f y , f 2 } is the velocity, = 
{likx.vky^kz} is the quasimomentum vector, and fi^g = £Lq{huj x + + &z) is the 
frequency (chemical potential in BEC arrays). 

By applying the quasicontinuum approximation 0, Q] we obtain for the first (F) 
harmonic /i = z/ = £ = lof the expansion ([3]) an equation which reads as 

a x d 2 x XF + a y dyXF + a z d 2 z XF ~ "oXf + 3U\xf\ 2 Xf = 0, (6) 
where Xf = Xi,i,i, and 

a = 2[J x cos(k x ) + Jycos(ky) + J z cos(k z )](l - Q ), 
otx = -J x cos(k x ), 

OLy = —Jy COS(ky), 

a z = -J z cos(k z ). (7) 

In Eq. (jSJ) the coordinates {x, y, z] = Si,i,i, where is given in Eq. (TJJ. In Eq. ([7]) the 

constants flo, k x , k y and fc 2 can be chosen as independent parameters. 

Since the coordinates {x,y, z] are mutually independent, Eq. (jSJ) can be integrated. So, 
finally, approximate soliton solutions of Eq. ([TJ) read as 



[^ )1 26 S ech3(-— )+^ ! _ 1 tanh3(-— )], 

m=x,y,z 



(8) 



where ^,±1 is the Kronecker delta with /3 W = siga(j w ao/a w ) and w = x,y,z. In Eq. (JHJ) 
the soliton widths read as 



Lbw = \fl w ao/a w (9) 

and 

L Dw = yJ-lwOLo/{2a w ), (10) 

where the relation j x + 7 y + 7^ = 1 should be satisfied. Besides 

©ill = k x m + k y n + k z p — fi^ijt (11) 

^1,1,1 = — 2^o [Jx cos(k x ) + JyCos(k y ) + J z cos(k z )}, (12) 

and 

vi,i,i = 2{J X sm(k x ), J y sm(k y ), J z sin(A; 2 )}. (13) 

In Fig. [T| it is shown an example of the motion of solitons following from initial conditions 
given by Eq. (jSJ). 

We note that other soliton solutions with embedded vorticity can be calculated by using 
spherical coordinates (r,9,(f)). In that case Eq. ([6]) can be written as 

,2 D n D |2 1\ , a (Jl . 



(3r 2 R{\R\ 2 -l) + d r (r 2 R) + 



d ( ((l - C)d ( R) + dlR/(l -C 2 ) = 0, (14) 



where 



R= v WJ^~oXf, (15) 

r 2 = a (x 2 /a x + y 2 /a y + z 2 /a z ), (16) 

C = cos(0), (17) 

= tan~ 1 [ v /7i~r/ (y/a^y)], (18) 

<f> = cos [y/aoz/(y/a^r)]. (19) 



By using the ansatz R = F(r)YJf(9, </>), the equation for R can be reduced to 



Tr ( r2 t f ) + Pr2fif2 ~ 1} " V{V + 1)f = °" (20) 
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FIG. 3. (Color online) Superposition of snapshots similar as in Fig. [IJ but with k = 0.2. Arrows 
in yellow show the path and direction of motion. At t = [labels A, B, and C (red color)] the 
initial conditions follow from Eq. (|2ip . where (A) k = {tt, tt, — 1.05^}, f^o = 1.05, v = /i = 
1 [extended dark vortex plotted in Fig. [2jb) (orange double-dot-dash line) truncated with the 
function 0.5(1 - tanh(r - 4)]], (B) Q = 1.01, v = \i = 0, and (C) Q = 1.05, v = n = 0. At t = 10 
[label D (green color)] the collision can be observed. And at t = 24 [labels A', B', and C (blue 
color)] the vortex and solitons after the collision can be observed. Other parameters are given in 

Fig. m 

Here, / = \Y£ \F(r) and YJf is a spherical harmonic function of degree v [y > 0) and order 
M (M < v)- in Eq. (l20l) the parameter (3 — 1 (/3 = —1) for ao« w > (a a w < 0) with 
w = x,y, z. 

Finally we obtain the solution 

^ = ymmoM m exp(z01 ' 1 - l) ' (21) 

where the embedded vorticity is given by the order parameter fi. In Eq. ( 12T1) the behavior 
of the function / is governed by Eq. (120|) . Notice that for v = a solution of Eq. (1201) is 
/ = 1. For given ^ and /3 values in Eq. fl20|) . other unbounded and bounded solutions can 
be numerically calculated with a shooting method. Some of these / solutions are plotted 
in Fig. Efb). Moreover, an example of the motion of ip m ,n,p functions given in Eq. ( |2TT) is 
shown in Fig. [31 
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III. RESULTS 



The split-step Fourier method has been used to solve numerically Eq. ([I]) for a lattice 
size max(m x n x p) = 64 x 64 x 64 with periodic boundary conditions. 

In Figs. [T] and [3] superpositions of snapshots at different time values of a soliton-soliton- 
soliton and a vortex-soliton-soliton collision are presented. The plots consist in p m ,n,p isosur- 
faces (surface for which p m)n . p is constant). The initial conditions in Figs. [Hand [3] [labels A, 
B, and C (color red)] follow from Eqs. dHJ and (I2TI) . respectively. The direction of motion 
of solitons and vortices is given by the the normalized velocity vector, 

d = vi >1| i/|vi, 1) i|, (22) 

where v^i i is given by Eq. ffTB"]) . 

In Figs. Q] and [3J we observe that solitons and vortices when propagating undergo a self- 
defocusing instability. This inestability follows from the interplay between nonlinearity and 
discreteness. Notice that nonlinearity tends to localize the wave while the discrete diffraction 
tends to spread it out. For moving solitons the self-defocusing instability corresponds to 
the case when diffraction effect is stronger than the nonlinear localization effect. This self- 
defocusing process manifest itself as a soliton-shape broadening accompanied by exponential- 
like amplitude decay. An example of this amplitude decay can be observed in Fig. [2](a), 
where the evolution of the absolute p m ,n,p maximum (p m ax) of the collision in Fig. [I] is 
plotted. 

Notice that in Fig. Hl^a) Pmax peaks when the collision in Fig. [1] [label D (green color)] 
occurs. After collision the solitons further broaden in shape, as can be observed in Fig. 
CD [labels A', B', and C (blue color)]. Besides, for comparison in Fig. [2](a) the amplitude 
Pmax of the individual solitons in the absence of collision is plotted. We can observe that 
before and after collision the pmax value coincides with that of individual solitons. Moreover, 
we can observe that the peak amplitude during collision is much higher than three times 
the amplitude of the individual solitons. It is because during collision the nonlinear effect 
becomes strong, so a transitory self-focusing process takes place. 

The self-defocusing process in solitons depends on their form and direction of motion, 
which in turn are governed by k. As in two dimensional systems (si, 9], solitons moving 



S 



along the diagonal directions, 



d 



< 



' {±1,±1,±1}/V3 
{±1,±1,0}/V2 
{±1,0,±1}/V2 
k {0,±1,±1}/V2 



(23) 



are less prone to self-defocusing than those moving along the main-axes directions, 



Notice that it is possible to control the defocusing process so that solitons moving in different 
directions and velocities can have the same decay rate. For example, solitons in Fig. [I] have 
similar decay rates, as shown in Fig. Ufa). This can be easily achieved by first choosing 
with the help of k the wanted directions and velocities of motion and then with the help of 
Qo, included in ao in Eqs. (JZj), the same initial amplitudes. Here we use the fact that the 
amplitude of a soliton is proportional to its decay rate 0, Q] . Notice that for similar decay 
rates we obtain solitons with different sizes and velocities, as shown Figs. [I] and [3j 

In Fig. [3] we plot the initial soliton and vortex solutions [labels A, B, and C (color red)] 
following from Eq. (|2~T]) . In particular, the motion of a vortex shell (A — > D — > A' in Fig. 
[3]) along its azimuthal axis (p axis) with vorticity /i = 1 is shown. The soliton dynamics 
is similar as in Fig. (TJ however small perturbations of the soliton shape are observed after 
the collision (labels B', and C in Fig. [3]). The broadening and breaking-apart of the vortex 
in Fig. [3] is due to the modulation instability and can be also observed in the absence of 
collisions. This instability effect on vortices is more pronounced for higher values of the 
vorticity (/i > 2). 

We remark that the center of mass of both solitons and vortices move with a constant 
velocity whose magnitude measured in simulations is very well predicted by Eq. ( 1T31) . On 
the other hand, though not seen in Figs. [1] and [31 the presence of low-intensity radiation 
tails can be detected for both soliton and moving vortices. 

A question that emerges from the results above is how to mitigate the self-defocusing 
effect. In order to tackle this problem we investigate the effect of a magnification in the 



{±1,0,0} 



d= <^ {0,±1,0} 



(24) 



{0,0, ±1} 
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form of a 3D thin-lens phase, i.e. 

ipm,n, P -> ipm,n, P exp(-ir 2 /(4:T)). (25) 

In the present analysis ipm,n,p solutions follow from Eq. (1211) . where the radial coordinate r 
is given by Eq. f|T6|) . In Eq. fl25l) . due to the discreteness of the system, the parameter T 
does not correspond exactly to the focal length, as it happens in continuous systems. 

In Fig. 11(a) we consider three solitons moving in the positive direction of the m axis, and 
their individual p max values are plotted in Fig. BJb). In Fig. H(a) the transversal separation 
distance between the solitons is large enough to avoid soliton-soliton interactions. The initial 
soliton forms, ip [Labels C, A, B (red color) in Fig. Ufa)], are identical to each other but 
distinct from the T value in the imposed thin- lens phase. In order to compare the soliton 
behavior, the level value of the isosurfaces in Fig. 11(a) was defined to be fixed to some 
initial value, as shown in|H(b) (red line). So, solitons with p max value above this level appear 
plotted in Fig. H(a). Notice, e.g., that in the soliton evolution A — > A' — > A" in Fig. @Ja) 
the soliton shape vanishes at A" because its p ma x value [black line in Fig. H(b)] decays below 
the level value [red line in Fig. IH(b)] of the isosurface for t > 12. 

Three different T values have been considered in Fig. H](a), namely T = oo, 0.2, and 
0.1. The value T = oo (A — > A' — > A") corresponds to the case where magnification is 
negligible. The value T = 0.2 (B — > B' — > B") corresponds to a case where for scale of time 
considered the soliton amplitude remains nearly constant [see Fig. 11(b), blue dashed line]. 
However, it is to remark that for larger time scales a self-defocusing process is observed and 
unavoidable. The value T = 0.1 (C — > C' — > C") corresponds to a case where a focusing due 
to the thin- lens phase is immediately followed by a self-focusing effect, as can be observed 
in Fig. S](b) (purple dotted line). Besides, a strong radiation tail can be also observed in 
Fig. ffl[a) (C"). 

For further comparison in Fig. 11(b) the case of T = 0.05 (green dot-dashed line) is also 
plotted. This case is similar to the case T = 0.1 where the first peak corresponds to a 
focusing due to the imposed thin-lens phase and the second peak is due to a self-focusing 
effect. Notice that after the second peak the soliton amplitude strongly decays and a long 
radiation tail, more in the fashion of Fig. H(a) (C"), can be also observed. 
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FIG. 4. (Color online) a: Superposition of snapshots of p m ,n,p isosurfaces at different time values 
of 3 solitons with different thin-lens-phase factors. Here, the isosurfaces are defined as the fix value 
Kmax„i W Pm,n,p(t = 0) and k = 0.5. Arrows (in yellow) show the path and direction of motion. 
At t = [labels A, B, and C (red color)] the initial conditions follow from Eq. (JSj) , where (A) 
T = oo, (B) T = 0.2, (C) T = 0.1. At t = 12 [labels A', B', and C (green color)] and t = 24 
[labels A", B", and C" (blue color)] solitons move with different maximums. Other parameters: 
k = 0.95f {1,0,0}, n = 1.001, and v — fi — 0. b: Evolution of Pm.n,p maximum of individual 
solitons. T = oo (Black solid line), T = 0.2 (blue dashed line), T = 0.1 (purple dotted line), and 
T = 0.05 (green dot-dashed line). The straight red line corresponds to the isosurface level of the 
plots in panel (a). 
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IV. CONCLUSIONS 



We have studied, for the first time, the motion of solitons and vortices in the 3D-DNLSE. 
In the tight binding limit this is the most simple model for studying the effect of nonlin- 
earities in 3D crystal-like structures, such as 3D photonic crystals, metamaterials and 3D 
BEC arrays. The analytical results, which are supported by simulations, suggest that these 
moving excitations can appear or be excited for any finite value of the model parameters. 
This implies that solitary waves may play an important role in the design of 3D crystal-like 
structures, such as cloaking devices. On the other hand, we have observed that solitons 
survive collisions and their modulation instability can be delayed. This supports the idea 
that photonic crystals may be used as a 3D compact routing [l| for optical pulses. From 
the practical standpoint, it remains challenging the exact generation of those solitons stud- 
ied here. However, it does not preclude the fact that moving excitations observed in those 
systems can be analyzed as a superposition of solitons and/or vortices. 



[1] D. Shir, E. C. Nelson, Y. C. Chen, A. Brzezinski, H. Liao, P. V. Braun, P. Wiltzius, K. H. A. 

Bogart, and J. A. Rogers, App. Phys. Lett. 94, 011101 (2009). 
[2] M. S. Rill et al, Nature Materials 7, 543 (2008) . 

[3] Y. Lai, H. Chen, Z.-Q. Zhang, and C. T. Chan, Phys. Rev. Lett. 102, 093901 (2009). 
[4] A. Greenleaf, Y. Kurylev, M. Lassas, and G. Uhlmann, Phys. Rev. Lett. 101, 220404 (2008). 
[5] I. V. Shadrivov, A. A. Zharov, N. A. Zharova and Y. S. Kivshar, Photonics and Nanostructues 
4, 69 (2006). 

[6] P. Chak, J. E. Sipe, S. Pereira, Opt. Lett. 28, 1966 (2003) 

[7] A. V. Yulin, A. R. Champneys, and D. V. Skryabin, Phys. Rev. A 78, 011804(R) (2008). 
[8] H. Susanto, P. G. Kevrekidis, R. Carretero-Gonzalez, B. A. Malomed, and D. J. Frantzeskakis, 

Phys. Rev. Lett. 99, 214103 (2007). 
[9] E. Arevalo, Phys. Rev. Lett. 102, 224102 (2009). 
[10] P.G. Kevrekidis, B. A. Malomed, D. J. Frantzeskakis, and R. Carretero-Gonzalez, Phys. Rev. 

Lett. 93, 080403 (2004). 



12 



[11] R. Carretero-Gonzalez, P. G. Kevrekidis, B. A. Malomed, and D. J. Frantzeskakis, Phys. Rev. 
Lett. 94, 203901 (2005). 

[12] C. Chong, R. Carretero-Gonzalez, B. A. Malomed, and P. G. Kevrekidis, Physica D 238, 126 
(2009). 

[13] A. Neuper, F. G. Mertens and N. Flytzanis Z. Phys. B 95, 397 (1994). 
[14] E. Arevalo, Phys. Lett. A 373, 3541 (2009). 



13 



